Interactive Session: Exploring Urban Greenhouse Gas Emissions

About the Data

The Greenhouse gas And Air Pollutants Emissions System (GRA²PES) dataset at the GHG Center is an aggregated, regridded, monthly high-resolution (0.036 x 0.036°) data product with emissions of both greenhouse gases and air pollutants developed in a consistent framework. The dataset contains emissions over the contiguous United States covering major anthropogenic sectors, including energy, industrial fuel combustion and processes, commercial and residential combustion, oil and gas production, on-road and off-road transportation, etc. Carbon dioxide (CO₂) emissions are developed along with those of air pollutants including CO, NOₓ, SO₂, and PM2.5 with consistency in spatial and temporal distributions. Emissions by sectors are reported as column totals in units of metric tons per km² per month. Spatial-temporal surrogates are developed to distribute CO₂ emissions to grid cells to keep consistency between greenhouse gases and air quality species. The current version of GRA²PES is for 2021. Long-term emissions and more greenhouse gas species (e.g., methane) are under development and will be added in the future.

Requirements

  • Set up Python Environment - See setup_instructions.md in the /setup/ folder

Learning Objectives

  • How to use U.S. GHG Center STAC Catalog to access Greenhouse gas And Air Pollutants Emissions System (GRA²PES) data
  • How to visualize datasets using folium and perform zonal statistics
  • How to plot time series for Greenhouse gas And Air Pollutants Emissions System (GRA²PES) and analyze the results for a region of interest

Approach

  1. Identify available dates and temporal frequency of observations for the given collection using the GHGC API /stacendpoint. The collection processed in this notebook is the Wetland Methane Emissions, LPJ-wsl Model data product
  2. Pass the STAC item into the raster API /stac/tilejson.json endpoint
  3. Define the spatial region of interest
  4. Using plugins from folium to visualize two tiles (side-by-side), allowing time point comparison
  5. After the visualization, perform zonal statistics for a given polygon
  6. Plot monthly time series for GRA²PES

Setup

Import the required Python libraries by running the next cell.

from pystac_client import Client
import matplotlib.pyplot as plt
import matplotlib
import pandas as pd
import numpy as np
import hvplot.xarray
from hvplot import hvPlot
import xarray as xr
import rioxarray as rxr
import stackstac
import h5py 
import os
os.environ['USE_PYGEOS'] = '0'
import datetime
import geopandas as gpd
import folium
import folium.plugins
from folium.plugins import MousePosition
import cartopy.crs as ccrs
import cartopy.feature as cf
import cartopy.io.img_tiles as cimgt
import requests
import surfaceAreaGrid
from geopy import geocoders
# For reading COG to xarray
import stackstac
# For displaying image in a jupyter notebook
from IPython.display import Image, display
# For reading from S3 bucket
import boto3
pd.options.plotting.backend = 'holoviews'
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[15], line 23
     21 import cartopy.io.img_tiles as cimgt
     22 import requests
---> 23 import surfaceAreaGrid
     24 from geopy import geocoders
     25 # For reading COG to xarray

ModuleNotFoundError: No module named 'surfaceAreaGrid'
# Provide STAC and RASTER API endpoints
STAC_API_URL = "http://earth.gov/ghgcenter/api/stac"
RASTER_API_URL = "https://earth.gov/ghgcenter/api/raster"
# Fetching the collection from STAC collections using appropriate endpoint
# The 'requests' library allows a HTTP request possible
# Open catalog, print info and available collections
catalog = Client.open(STAC_API_URL)

# print catalog and collections titles
print('>> Title: ',catalog.title)
print('>> Collections:')
for collection in catalog.get_collections():
    print(collection.id)
>> Title:  US GHG Center STAC API
>> Collections:
ct-ch4-monthgrid-v2023
emit-ch4plume-v1
lpjwsl-wetlandch4-monthgrid-v1
vulcan-ffco2-yeargrid-v4
goes-ch4plume-v1
odiac-ffco2-monthgrid-v2023
tm54dvar-ch4flux-monthgrid-v1
lpjeosim-wetlandch4-monthgrid-v2
gosat-based-ch4budget-yeargrid-v2
lpjeosim-wetlandch4-daygrid-v2
gra2pes-ghg-monthgrid-v1
gosat-based-ch4budget-yeargrid-v1
micasa-carbonflux-daygrid-v1
sedac-popdensity-yeargrid5yr-v4.11
vulcan-ffco2-elc-res-yeargrid-v4
tm54dvar-ch4flux-mask-monthgrid-v1
casagfed-carbonflux-monthgrid-v3
micasa-carbonflux-monthgrid-v1
lpjeosim-wetlandch4-daygrid-v1
lpjeosim-wetlandch4-monthgrid-v1
oco2geos-co2-daygrid-v10r
epa-ch4emission-yeargrid-v2express
lpjwsl-wetlandch4-daygrid-v1
oco2-mip-meanco2budget-yeargrid-v1
oco2-mip-co2budget-yeargrid-v1
eccodarwin-co2flux-monthgrid-v5
odiac-ffco2-monthgrid-v2022
# Select GRA2PES collection: browse item_assets, spatial and temporal extent
collection = catalog.get_collection('gra2pes-ghg-monthgrid-v1')
collection
# Get items from this collection, examine temporal extent of each
items = collection.get_items()
for item in items:
    print(item)
<Item id=gra2pes-ghg-monthgrid-v1-202112>
<Item id=gra2pes-ghg-monthgrid-v1-202111>
<Item id=gra2pes-ghg-monthgrid-v1-202110>
<Item id=gra2pes-ghg-monthgrid-v1-202109>
<Item id=gra2pes-ghg-monthgrid-v1-202108>
<Item id=gra2pes-ghg-monthgrid-v1-202107>
<Item id=gra2pes-ghg-monthgrid-v1-202106>
<Item id=gra2pes-ghg-monthgrid-v1-202105>
<Item id=gra2pes-ghg-monthgrid-v1-202104>
<Item id=gra2pes-ghg-monthgrid-v1-202103>
<Item id=gra2pes-ghg-monthgrid-v1-202102>
<Item id=gra2pes-ghg-monthgrid-v1-202101>
# Explore one item (one month!): Check out properties, assets
item = collection.get_item('gra2pes-ghg-monthgrid-v1-202112')
item

Science topic: Compare urban and rural CO2 season cycles

Let’s compare an urban area to a rural area - for example, the Chicagoland area vs. central Illinois.

Read data from STAC into xarray Dataset

search = catalog.search(
    collections='gra2pes-ghg-monthgrid-v1',
    datetime=['2021-01-01T00:00:00Z','2021-12-31T00:00:00Z']
)
# Take a look at the items we found
for item in search.item_collection():
    print(item)
<Item id=gra2pes-ghg-monthgrid-v1-202112>
<Item id=gra2pes-ghg-monthgrid-v1-202111>
<Item id=gra2pes-ghg-monthgrid-v1-202110>
<Item id=gra2pes-ghg-monthgrid-v1-202109>
<Item id=gra2pes-ghg-monthgrid-v1-202108>
<Item id=gra2pes-ghg-monthgrid-v1-202107>
<Item id=gra2pes-ghg-monthgrid-v1-202106>
<Item id=gra2pes-ghg-monthgrid-v1-202105>
<Item id=gra2pes-ghg-monthgrid-v1-202104>
<Item id=gra2pes-ghg-monthgrid-v1-202103>
<Item id=gra2pes-ghg-monthgrid-v1-202102>
<Item id=gra2pes-ghg-monthgrid-v1-202101>
# Read data into an xarray Dataset
ds = stackstac.stack(search.item_collection(),epsg=4326).squeeze()
ds.sel(band='co2')[9,:,:].hvplot(x='x',y='y',cmap='Spectral_r',coastline=True,clim=(0,2000))
/srv/conda/envs/notebook/lib/python3.12/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
# Take a look at the dimensions of the data
print(ds.sel(band='co2').shape)
print(ds.sel(band='co2').dims)
(12, 947, 2188)
('time', 'y', 'x')

Select areas of interest using GEOJSON shapes and rioxarray

Read in AOIs

GEOJSONs generated via geojson.io

data_dir = f"{os.getenv('HOME')}/shared/data/use-case-2"   
central_IL = gpd.read_file(f'{data_dir}./atlanta.geojson').to_crs(epsg=4326)
chicagoland = gpd.read_file(f'{data_dir}./chicagoland.geojson').to_crs(epsg=4326)
import cartopy.feature as cf
cmap='Spectral_r'    # colormap definition
vmin=0               # minimum value on colorbar
vmax=1000            # maximum value on colorbar
fig1=plt.figure(figsize=(8,5))
ax1 = fig1.add_subplot(projection=ccrs.PlateCarree())
p1 = ds.sel(band='co2')[9,:,:].plot(vmin=vmin,vmax=vmax,cmap=cmap,ax=ax1)
# Rural is green
rural_color = '#00A676'
# Urban is pink
urban_color = '#eb7bc0'
central_IL.plot(ax=ax1,edgecolor=rural_color,linewidth=1.5,facecolor='none')
chicagoland.plot(ax=ax1,edgecolor=urban_color, linewidth=1.5,facecolor='none')
ax1.set_xlim(-95.4,-78.5)
ax1.set_ylim(32.6,43.04)
ax1.add_feature(cf.STATES,linewidth=0.5)

Clip CO2 field to AOIs using rioxarray

rio.clip has a few options for selection of grid cells relative to the specified geometry. Here we will use the default functionality. See documentation for more info.

rural_clip = ds.rio.clip(geometries=central_IL.geometry.values)
urban_clip = ds.rio.clip(geometries=chicagoland.geometry.values)
# We also need the latitude values reversed to utilize a function correctly in the next step.
rural_clip = rural_clip.reindex(y=rural_clip.y[::-1])
urban_clip = urban_clip.reindex(y=urban_clip.y[::-1])
import cartopy.feature as cf
fig1=plt.figure(figsize=(8,5))
ax1 = fig1.add_subplot(projection=ccrs.PlateCarree())
p1 = rural_clip.sel(band='co2')[9,:,:].plot(vmin=vmin,vmax=vmax,cmap=cmap,ax=ax1)
p2 = urban_clip.sel(band='co2')[9,:,:].plot(vmin=vmin,vmax=vmax,cmap=cmap,ax=ax1)
# Rural is green
rural_color = '#00A676'
# Urban is pink
urban_color = '#eb7bc0'
central_IL.plot(ax=ax1,edgecolor=rural_color,linewidth=1,facecolor='none')
chicagoland.plot(ax=ax1,edgecolor=urban_color, linewidth=1,facecolor='none')
ax1.set_xlim(-95.4,-78.5)
ax1.set_ylim(32.6,43.04)
ax1.add_feature(cf.STATES,linewidth=0.5)

Calculate and apply area weights for lat/lon grid cells

Emissions are in tons km-2 month-1, but our grid cells are in lat/lon space - calculate how many km are in each grid cell, so we can calculate tons month-1 for each AOI.

# This function returns m2, and we need km2. Note the /1000./1000.
rural_weights = surfaceAreaGrid.surfaceAreaGridd(lat_centers=rural_clip.y.values, \
    lon_centers=rural_clip.x.values,ret_area=True)/1000./1000.
urban_weights = surfaceAreaGrid.surfaceAreaGridd(lat_centers=urban_clip.y.values,\
    lon_centers=urban_clip.x.values,ret_area=True)/1000./1000.   
# Apply these weights - simple as a multiplication
rural_weighted = rural_clip.sel(band='co2')*rural_weights
urban_weighted = urban_clip.sel(band='co2')*urban_weights
# Take the mean of emissions over AOI grid cells
rural_means = rural_weighted.mean(dim=['x','y'])
urban_means = urban_weighted.mean(dim=['x','y'])
# Reformat some times for plotting purposes
times = [sd[0:10] for sd in rural_means.start_datetime.values]
times = [datetime.datetime.strptime(sd,'%Y-%m-%d').strftime('%b') for sd in times]

Plot urban vs. rural CO2 emissions time series for 2021

fig,ax1 = plt.subplots()
# Rural first
ax1.plot(times[::-1],rural_means.values[::-1],marker='o',color=rural_color,label='Chicago',lw=2.5)
#ax1.set_ylabel('Rural',color=rural_color,fontweight='bold')
# Urban next
#ax2 = ax1.twinx()
ax1.plot(times[::-1],urban_means.values[::-1],marker='o',color=urban_color,label='Atlanta',lw=2.5)
#ax2.set_ylabel('Urban',color=urban_color,fontweight='bold')
plt.legend()
plt.title('Mean urban CO$_2$ emissions, 2021')